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ABSTRACT 

Given  a  non-singular  matrix  then  factorlze  It  Into  the  product 
of  two  matrices.   Form  a  new  matrix  by  multiplying  the  two  factors 
together  in  reverse  order.   Repeat  the  process  and  so  produce  a 
sequence  of  matrices  all  similar  to  the  given  matrix.   This  is  the 
essense  of  an  LR  method  and  for  several  factorization  techniques 
the  sequence  converges,  quite  generally,  to  a  triangular  matrix 
from  which  the  eigenvalues  may  be  read  off. 

The  LR  transformations  discussed  here  preserve  the  form  of  a 
Hessenberg  matrix  (  a.  .  =  0   if   i  >  j+1  )  and,  for  such  matrices, 

the  work  involved  in  the  transformation  from  one   NX  N  matrix  to 

2  ^ 

the  next  is  proportional  to   N   as  against   N^  for  full  matrices. 

Using  this  and  other  devices  the  latest  versions  of  the  transform- 
ation are  among  the  best  available  eigenvalue  methods  for  general 
matrices . 

The  evolution  of  these  algorithms  can  be  traced  back  to 
classical  function  theory,  each  method  being  more  feasible  numer- 
ically than  its  predecessor. 
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THE  DEVELOPMENT  AND  USE  OF  METHODS  OF  LR  TYPE 

by 
Beresford  Parlett 

1 .   Introduction 

The  rapid  proliferation  of  numerical  methods  in  linear  algebra 
has  made  the  family  tree  of  their  development  almost  unreadable. 
Therefore  it  Is  an  unexpected  pleasure  to  find  a  section  of  the 
tree  whose  evolution  can  be  viewed,  without  too  much  distortion, 
as  a  simple  sequential  process,  each  development  removing  a  defect 
of  the  previous  method.   This  paper  is  written  from  the  point  of 
view  of  someone  who  seeks  a  "good"  computational  way  of  finding 
all  the  eigenvalues  of  a  general  matrix.   This  viewpoint  has 
shaped  the  whole  article  and  caused  the  exclusion  of  interesting 
but  irrelevant  facets  of  some  of  the  topics. 

The  aim  of  the  paper  is  to  present  the  ideas  and  character- 
istics of  the  various  algorithms  in  relation  to  each  other  so 
that  the  subject  can  be  seen  as  a  whole.   Ample  references  are 
supplied  to  sources  which  give  full  details. 

One  way  to  find  the  eigenvalues  of  a  matrix  is  to  generate 
from  it  a  sequence  of  similar  matrices  whose  limit  is  triangular. 
The  methods  considered  here  are  of  this  type.   Each  term  is  factor- 
ized  into  the  product  of  two  matrices.   The  next  term  is  then  the 
product  of  these  factors  in  reverse  order.   In  addition  these 
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methods  can  be  regarded  as  developments  of  the  QD  algorithm  which 
Itself  can  be  seen  as  a  practical  way  of  using  Hadamard's  character- 
ization of  the  poles  of  a  rational  function  in  terms  of  the  coef- 
ficients of  its  Taylor  series  about  a  regular  point. 

The  sections  on  QD  are  based  on  Henrici's  interesting  article 
[10],  though  the  viewpoint  is  different  and  motivation  is  given  for 
the  introduction  of  certain  essential  quantities.   Rutishauser ' s  LR 
transformation  [l4]  is  suggested  by  an  aspect  of  his  QD  algorithm. 
The  QR  transformation  of  Francis  [6]  avoids  the  possible  "instability" 
of  the  original  LR  transformation  but  is  slower.   Wilkinson  uses  the 
LR  method  with  row  and  column  interchanges  to  produce  an  algorithm 
which  is  stable  and  requires  no  more  arithmetic  than  the  basic  LR. 
In  practice  this  variation  seems  both  fast  and  accurate  but  no  con- 
vergence proof  has  been  published  yet  nor  is  it  apparent  how  to 
adapt  the  technique  to  give,  economically,  the  complex  eigenvalues  of 
real  matrices . 

The  final  section  describes  recent  work  of  Ortega  and  others 
to  produce  a  rival  to  the  Sturm  sequence  method  of  finding  the 
eigenvalues  of  real  symmetric  tridiagonal  matrices. 
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2 .   Classical  Theory 

In  order  to  find  the  eigenvalues  of  a  matrix  A   it  is  com- 
pletely satisfactory  to  replace  A   by  A  +  si   (I   is  the  iden- 
tity) where   s   is  any  convenient  number.   Hence  there  is  no  loss 
of  generality  in  considering  only  non-singular  matrices.   Unless 
otherwise  stated  A   is  to  be  taken  as  a  complex   NX  N  matrix  with 
eigenvalues   A.   ordered  by  the  relation 

|Aj  >  U^l  >  ...  >  lA^I  >  0  . 

The  starting  point  for  the  investigations  presented  here  is 
the  casting  of  the  eigenvalue  problem  in  a  function  theoretic  form 
which  is  going  to  prove  most  fruitful.   Let  x*   denote  the  con- 
jugate transpose  of  the   N  dimensional  vector  x  .   For  (almost) 
arbitrary  vectors   x  and  y  we  choose  to  consider  a  certain 
rational  function  of  a  complex  variable   z  ,  regular  at   0  and  oo  , 
which  is  defined  by 

f(z)   -  y*(I-zA)"^x   . 
The  eigenvalues  of  A  are  the  reciprocals  of  the  poles   z. 
of   f  .   Moreover  since   f   is  regular  at  the  origin  it  is  fully 
determined  by  its  Taylor  expansion  there  (which  is  convergent  for 

f(z)   =  ^  a^z^  ,   a^   =  y^A'^x  ,   A°   =  I  . 
0 
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In  theory  the   a   are  easily  calculated  from  A   and  so  may 
be  taken  as  known.   The  problem  is  now  a  special  case  of  the  problem 
considered  and  solved,  in  principle,  by  Hadamard  at  the  end  of  the 
nineteenth  century  [9]. 

Determine  the  singularities  of  a  meromorphic  function  from  the 
coefficients  of  its  Taylor  series  about  a  regular  point. 

In  l884  K8nig  proved  that,  if  the  only  singularity  on  the 
circle  of  convergence  of  f   is  a  simple  pole  at   z-,   then 


Z-,   as   n  — >  00  .   Hadamard  [9]  solved  the  general 


(a  /a  ^.  )  - 
n     n+1 

problem  in  I892  with  the  aid  of  certain  determinants,  now  usually 
called  Hankel  determinants,  defined  as  follows 


H- 


n 


xn 


1,   H   =  det 

'    k 


n 


n+1 


.  a 


a 


n+1 


a 


n+2 


n+k-1 


n+k 


^n+k-1  ^n+k 


a 


n+2k-2 


(n  =  0,1,2, 
k  =  1,2,.. 


The  basic  result  is  contained  in  the  following 

00 
LEMMA  1.   Let   f(z)  =  ^        a  z"^  =  f  (z)  +  r(z)  ,  where   f    is 

11-1  -1 

regular  analytic  in  \z\   ^  L  -^    and  r  is  rational  with  poles   A. 

satisfying 


^1  >  \\\    > 


>  U    >  L  >  0 
—  '  m' 


(with  the  convention  that  a  pole  of  order  p   is  regarded  as   p 
concurrent  simple  poles)  then  as   n  — >  00 


H' 


n 
"m 


C(A 


,n 


1 


J.J''  [1   +  o{(lA^)'^)}  , 


where   C   is  a  non-zero  constant  Independent  of  n  . 

This  is  proved  by  Henrici  in  [10]  for  the  case  in  which  all  the 
A.   are  distinct.   The  general  case  is  covered  by  Golomb  in  [8]. 

The  essence  of  the  proof  lies  in  the  fact  that   a   is  of  the 
form 


^n  -  ^r'l  +  •••  +  Vn 


a._   =  B,  AV  +  .  .  .  +  B^At;  +  B. 


and  so  to  find  A-,A^...A   it  is  only  necessary  to  find  some  function 

of  the   a  ,a  ,-,,,..   which  is  of  the  form 
n'  n+1 ' 

(A,  ...  A  )'^  F 
1      m 

where   F  may  depend  on  anything  but  n  .   It  turns  out  that   H   is 
Just  such  a  function. 

From  the  lemma  follows 

COROLLARY  1.   For  large  enough  n  ,   H^  ^  0  . 

COROLLARY   2.      As      n  ^-  oo  ,      (  H^"^^/H^  ) -^  A.  .  .  .  A      . 
—  \^  m     /    m  /  1  m 

So  classical  function  theory  gives  the  solution  that  the  k  th  eigen- 
value of  A   is  given  by 


n->m  \         /  ' 


n->co 

As  Henrici  points  out  this  is  unsatisfactory  from  the  n\imerical 
point  of  view  as  there  is  no  simple  way  to  compute   H.  .   Aitken 
overcame  this  defect  very  neatly  and  produced  a  feasible  numerical 
procedure.   Rutishauser,  with  the  QP  algorithm,  showed  how  to  by- 
pass the  computation  of  the  Hankel  determinants  altogether  with  a 
resulting  reduction  in  the  amount  of  computation. 
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3.   Altken's  Scheme 

During  the  1920' s  Altken  worked  on  the  problem  of  finding  all 
the  zeros  of  a  polynomial  by  generalising  Bernoulli's  method  [1] 

and  was  led  to  study  the  Hankel  determinants.   He  discovered  a 

n 
relation  which  allows  the  recursive  calculation  of  the  H,   from  the 

a.  .   This  is  most  easily  shown  by  arranging  the   H,   in  the  follow- 
ing table  : 


H? 

H^ 

H° 

H? 

4 

4 

H? 

4 

4 

K 

4 

4 

4 

4 

4 

(H  -  table) 


Note  that  H-,  =  a   and  so  the  first  two  columns  are  known. 
1    n 

In  the  case  under  consideration   f(z)   has  only  N  poles  and 

so  the   H-table  has  only  N+l   columns.   This  follows  from  a 

classical  result  (see  Bieberbach  [3],  p.  319). 

LEMMA  2.   Let   f  (z)  =  Y~   a  z^   =    (b  +.  .  .+b  zP)/(c  +.  .  .+c  z^ ) 
^ n      '  o      P      o      q 

be  a  rational  function  (c   ?^  0,  c   ^  0)  .   Then 
_ Q  /      q         

^q+l      =   0  »        n  >  max  (0,p+l-q)  . 

Conversely  let   p   be  an  integer  such  that  E^   ^  0    ,      H^^  .,  =  0 
q  ^        '        q+1 

for  n  >  p+l-q  .   Then  f   is  of  the  above  form. 
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In  the  present  case  the  first  half  of  the  lemma  is  just  a  con- 
sequence of  the  Cayley-Hamilton  theorem  which  yields 

M        N-1 
A^  =  -  ;^  d.A^ 
0   ^ 

where  the   d.   are  the  coefficients  of  the  characteristic  polynomial 

of  A  .   Hence  for  k  =  0,1,2,... 

1=0  1=0 

Thus  for  n  =  0,1,2,...   the  last  column  of  H^,-i   is  a  linear  com- 

Ir 

bination  of  the  previous  columns.   Therefore  H,,  .,  =  0  . 

The  ratios  of  successive  elements  In  the   k  th  column  of  the 
H-table  converge  to  the  product  of  the  largest   k  eigenvalues  by 
Corollary  2  of  Lemma  1,  provided  that   |  A,  |  >  I  A,  .,  |  .   By  Corollary 
1  the   H-table  ultimately  exists,  i.e.   H,  ?^  0  ,   k  j<  N  ,   n   large 
enough.   However  it  is  possible  for  an  early  term  to  vanish  and  this 
would  prevent  the  recursive  computation  of  the  table  from  previous 
terms.   This  possibility  will  be  Ignored  here. 

The  relation  between  neighboring  Hankel  determinants 

The  relation  which  Aitken  found  for  himself,  [2],  but  which  had 
been  known  to  Hadamard,  is 

/„nx2   „n-l„n+l   ^n-l„n+l     ^ 

^V    -  \  \    +  ^k+i\-i   =  °  • 

This  is  best  remembered  by  taking  any  H,   as   P  and  labelling  its 
nearest  neighbors  in  the  H-table  topographically  as  N,  S,  W,  and 
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E  .   With  such  a  stencil  the  relation  above  can  be  written  as 

P^   =  NS  -  WE  .  (R) 

This  formula  can  be  used  In  two  different  ways: 

(a)  Start  with  the  first  two  columns  and  work  from  left  to 
right  In  the   H-table,  adding  one  column  at  a  time  by  repeated  use 
of   E  =  (NS-P^)/W  . 

(b)  If  the  first  two  diagonals  are  known  one  can  compute  the 
lower  diagonals  successively  working  along  the  new  diagonal  from 
left  to  right,  using   S  =  (P^-WE)/N  . 

Now  (b)  has  the  difficulty  of  finding  the  first  two  diagonals 
and,  although  he  showed  how  to  do  this  for  a  polynomial  whose  coef- 
ficients are  given,  Altken  In  [2]  favoured  method  (a). 

For  the  computation  of  eigenvalues  the  generation  of  the   H- 

table  by  (a)  or  (b),  though  feasible,  still  leaves  a  lot  to  be 

n      n 
desired.   With  (a)  It  Is  necessary  to  compute   a   =  H,  =  y*A  x 

for  as  many  n  as  are  necessary  to  obtain  convergence  of  H,   /H 

to  the  desired  degree  of  accuracy.   With  (b)  there  Is  the  necessity 

of  computing  a.   for  1  =  0,l,...,2N-2  and  then  the  determinants 

H,   and  H,   for  k  =  1,2,...,N  .   If  Gaussian  elimination  were 
k        k  >  >    » 

used  then,  whilst  calculating  H^  and  H^  ,  all  the  other  H. 
and  H.   (1  =  1,2,...,N-1)  ,  being  leading  principal  minors,  could 
be  found.   This  technique,  however,  would  not  permit  the  use  of 
Interchanges  In  the  Gaussian  elimination  and  so  extra  precision 

■X- 

would  be  required  to  ensure  adequate  accuracy.    Thus  both  (a)  and 

In  all  2W     multiplication  for  the  a.   and  2/3  N-^  double 
precision  multiplications  for  the  H.   and  H.  . 
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(b)  demand  a  formidable  amount  of  computation  before  beginning  to 
find  the  remainder  of  the   H- table. 

It  was  Rutlshauser' s  achievement  to  change  this  scheme  Into  a 
related,  but  far  from  obvious,  form  which  Is  numerically  preferable 
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4.   The  QD  Algorithm 


What  changes  could  be  made  in  the   H-scheme  which  might  reduce 

the  amount  of  Initial  computation?   There  is  one  natural  observation 

n 
the   H,   themselves  are  of  less  interest  than  the 

„n+l  n 
n  ^  \     \-l 

k  k-1 

since,  as  stated  in  section  2,    for   |  A-,  |  >  |  Ap  |  >  ...  >  |A:^j|  ,  as 
n  — >  00 

^k  -^  \  • 

As  Henrici  says  in  [10]:   "It  is  remarkable  that  in  the  comput- 

n  n 

ation  of  the   q,  ,  the  determinants   H,   do  not  have  to  be  used  if 

a  set  of  auxiliary  quantities  is  Introduced."   Their  introduction  is 

presented  as  a  "fait  accompli"  by  Rutishauser  in  [15].   Certainly 

this  set  of  quantities  arises  naturally  in  looking  at  the  continued 

fraction  expansion  of   f(z)   (see  [10],  p.  36)  but  this  author  feels 

that  a  motivation  for  the  choice  of  the  auxiliary  quantities,  called 

n 
k 

The  key  is  to  recast  the  basic  recursion  relation  (R)  so  that 


n 
e,  ,  can  be  made  in  the  present  context 


it  shows  explicitly  the  connection  between   qf  and  qP"*"   .   To  do 
this  clearly  consider  a  typical  portion  of  the  H-table  obtained  by 
placing  the  following  stencil  on  the  table  so  that   D  is   on  B^    . 
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B 


A 

C 

E 


D 

F 
H 


G 
I 


Then 


n   _  


AF 
CD 


n+1 
Ik 


CH 
EF 


To  Include  all  these  elements  two  applications  of  (R)  are  needed, 
centred  at   C  and  F  respectively,  namely 


C' 


=     AE  -  BD  , 


DH  -  EG 


Now  turn  the  first  terms  on  the  right  hand  sides  Into  q,^  and  q. 


^k 


respectively,  obtaining 


CF 

AF 

BF 

CF 

CH 

CG 

DE   " 

'-      CD  ■ 

CE  ' 

DE   " 

'   EP 

DF 

(rr: 

Reference  to  the  stencil  shows  that   BF/CE  and   CG/DF  are  of  the 
same  form  and  naturally  ask  to  be  named.   If  CG/DF 


(=  ^k-i^S+i/^X""^  ^  ^^  ^^11^^ 


e|J   then  BF/CE  will  be   eJJ^]^ 


Transposing  the  second  terms  on  the  right  hand  sides  of  (RR)  gives 
the  following  simple  relation  between  the  q  and  the  e  , 


n  ,   n 
^k  +  ^k 


n+1  ,   n+1 
^k   ^  ^k-1 


(R^) 


Another  relation  between  the  q  and  the   e   is  required  to 
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determine  all  the   q,   and  e,   from  the   q.,  (=  a   -, /a  )  .   It 

^k        k  ^1     n+1^  n 

follows  directly  from  the  definitions  of  the  q  and  e  .  A  simple 
derivation  comes  from  observing  that 

n   n   _  AF  CG   _  AG     n+1 .  n+1      CH  BF   _   BH 

^k   k      CD 'DP   "   „2  '   ^k    ^k-1   ~   EP'CE   "   „2  ' 

L)  Ji 

Reference  to  the  stencil  shows  that  the  right  hand  terms  are  of  the 
same  form  but  located  in  different  diagonals,  namely  n  and  n+2  . 
By  changing   n  and   k   so  that  both  terms  become   Cl/F    (in 


diagonal   n+1  )  gives 


n+1   n+1      n    n  /  0  \ 

^k   -^k     =   ^k+l'^k   •  (^2^ 


The  QD  Scheme 

Formula  ( R-,  )  shows  that   e,   is  the  difference  of  the  quotients 
q,     and   q,   modified  by   e,  _-,  .   This  indicates  a  little  why 
Rutishauser  named  the  following  triangular  table  of  q's   and  e's 
the  Quotient-Difference  Scheme. 


0  .  e° 

o 

^1 

^1  o 


o 


^1  .  Q' 


0   =  e  e^  e_ 

o  1  2 


2£  1  o 

^1  ^^2  ^3 


o       5  2  -^     ">\  1 

0  =  e  e,  e^  e 


o  1\+     ^"2  "3 

3       ^\  2  ^         1  o 

^1  ^2  ^J>  ^4 
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Typical  cases  of  (R-,  )  and  (Rp)  are  shown  plctorlally  In  the  table 
and  explain  why  E.  Stlefel  named  them  the  rhombus  rules.   (R, )  yields 
rhombl  centred  on  q-columns,  ( Rp )  yields  rhombi  centred  on  e- 
columns . 

By  alternate  application  of  the  rules  the  table  may  be  gener- 
ated, in  theory,  from  either 

(a)  the  first   q-column  moving  right,  or 

(b)  the  first  diagonal  moving  down. 

Now,  however,  method  (a)  is  quite  impractical  since  rule  (R-,  )  shows 
that 

n      n+1    n 
®1   "  ^1    ~  ^1  ■ 

For   |A-,  I  >  |Ap|  ,   q?  — >  X.  ,   as   n  — >  oo  .   As   n   increases 
e^  becomes  the  difference  of  two  almost  equal  numbers  and,  for 
computation  with  a  fixed  number  of  digits  for  each  number,  it  will 
have  fewer  and  fewer  significant  figures.   Thus   (e-j^   /e^ )   and 
hence   qp  will  have  little  or  no  accuracy.   Briefly,  one  says  that 
method  (a)  is  "unstable."   Method  (b),  on  the  other  hand,  uses  (R-^) 
in  the  form 

n+1      n  ,  ^n     „n+l 
^k    =  ^k  ""  ^k+1  -  ^k-1 

which  is  ultimately  completely  safe  when  eJJ_^^  and  e^_^  —>   0  . 

Difficulties  in  calculating  the  QD  scheme 

The  danger  in  method  (b)  is  that  in  the  early  stages,  for  small 
n  ,  large  values  of  ^w'^^  ^^y  occur  which  will  cause  loss  of  signi- 
ficant digits  in  the   q!^"^"^  .   However  with  the  QD  scheme  there  is  no 
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alternative  but  to  use  (b),  the  so-called  progressive  method. 
Rutishauser  in  [13]  has  analysed  the  effect  of  one  small   q,     and 
the  resulting  large  value  of  e,     (in  using  (Rp)).   He  observed 
that  only  negative  powers  of  e,     occur  in  diagonal   n+3  and 
below.   Thus  the  disturbance  affects  only  the  next  two  diagonals 
and  he  proposed  some  auxiliary  modifications  to  be  made  in  calculat- 
ing some  of  the  elements  in  diagonals   n+1   and  n+2  .   For  diagonal 
n+3  and  below  the  straightforward  application  of  the  rhombus  rules 
is  again  sufficient.   These  precautions  however  are  difficult  to 
incorporate  into  an  automatic  program  which  must  allow  for  the 
occurrence  of  several  large   e   values  in  different  columns  of  suc- 
cessive diagonals. 

This  same  danger  will  appear  again  In  the  LR  transformation  and 
recent  developments  can  be  regarded  as  techniques  to  overcome  this 
difficulty.   See  also  the  discussion  by  Wilkinson  [l8]. 

Although  as   n  — >  oo  ,   q,  7^->  0  by  assiomptlon,  it  is  still 
possible  to  have  q,  =  0  for  some  small  (finite)  value  of  n  .   If 
this  happens  the  QD  scheme  cannot  be  completed.   Formally  the  scheme 
is  said  not  to  exist.   This  possibility  will  be  ignored. 

Calculation  of  the  first  diagonal 

At  this  point  it  is  worth  observing  that 

n+1 /„n      ^An+1 


^n     Tj-n+l /„n      *An+l  /  „.n 
^1   "    1   /  1   "  y^   x/y*A  X 


and  so  the  descent  of  column  1  of  the  QD  scheme  corresponds  to  the 
well  known  power  method  [17]. 
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How  can  the  first  diagonal  of  the  QD  scheme  be  found  indirectly? 
The  answer  is,  surprisingly,  by  use  of  Lanczos'  method  of  "minimized 
iterations",  [11],  with  initial  vectors   x  and  y   (from  the  defin- 
ition of   f(z)   in  section  2).   Thus,  as  Henrici  points  out,  the  QD 
scheme  links  the  power  method  to  the  method  of  Lanczos. 

The  reader  is  referred  to  [11]  and  [10]  for  proofs  and  details 
of  the  algorithm.    It  produces  a  tridiagonal  matrix  similar  to  A 
and  of  the  form  (for  N  =  5  ) 


a. 


a^  P^ 


"il  Pi| 


a. 


By  the  LDU  theorem  ([4],  p.  20),  J  may  be  written  as  the  product  LR 
where 

'"1   e^ 


^1 

1 


q^ 


^3 

1   q 


4 


R 


"2 

1   e 


3 

1   e 


4 


It  is  easily  verified  that  for  k  =  1,2,...,N  , 


°'k  =  ^k  +  ^k-1  ' 


P, 


qi  ©1 

^k  k 


e^  =  0  , 


^N  =  ° 
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Henrlcl  shows  ([10],  p.  33)  that  these  q,   and  e,   are  the  elements  of 
the  first  diagonal  of  the  QD  scheme.   In  fact  he  shows  that  the 
n  th  diagonal  can  be  obtained  In  the  same  way  by  using  Lanczos' 
method  with  initial  vectors  A  x  and  y  .   In  particular 


a-,   =  y*A   x/y*A  x   =   q-. 


n 
'l 


as  it  should, 

Wilkinson  [l8]  and  Strachey  and  Francis  [l6]  have  examined  the 

numerical  reduction  of  a  full  matrix  A   to  tridiagonal  form  J  . 

The  recommended  method  is  to  reduce  A   to  an  intermediate 

(Hessenberg)  form  H  ,  indicated  here  for  N  =  5  , 

X  X   0   0   0 

X  X  X   0   0 

H   =    X  X  X  X   0 

X   X   X   X   X 
X   X   X   X   X 


This  reduction  may  be  performed  safely  and  accurately  in  single  pre- 
cision arithmetic  by  elimination  methods  if  row  and  column  inter- 
changes are  used.   It  requires   5/6  N-^  multiplications.   The  reduc- 
tion of  H  to   J  requires   1/6  N-^  multiplications  but  is  "un- 
stable" and  should  be  done  in  double  precision.   Comparison  with 
section  3  shows  that  the  first  diagonal  of  the  QD  scheme  can  be 
obtained  with  approximately  1/3  the  computation  required  for  the 

first  two  diagonals  of  the  H- table. 

I   I    I   I         II         n  n 

When   A-,   >   A„   >  ...  >   A, J   then  q,  — >  A,  ,   e,  — >  0  as 
'1'    '2'  'N'         ^k     k'k 


n 


GO 


When  there  are  eigenvalues  with  the  same  modulus  then 
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the  corresponding  q  and  e   columns  may  not  converge.   Fortunately 
this  does  not  Impair  the  usefulness  of  the  algorithm.   Discussion  of 
what  to  do  in  this  case  Is  most  simply  presented  In  connection  with 
the  LR  transformation  and  will  appear  In  the  next  section. 
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5.   The  LR  Transformation 


In  section  4  It  was  shown  that  the  elements   q.   and  e.   of 
the  first  diagonal  of  the  QJ)  scheme  are  obtained  from  the  LR  decom- 


position of  a  certain  trldiagonal  matrix  J  ,    now  to  be  called  J   , 


whose  non-zero  elements  In  row  1  are   1,  a.,  P> . 


This  Implies 


that,  for  k  =  1,2, 


,N 


o  ,   o 


k-1 
o  o 


a 


=  P 


o 

^° 

N 


0  , 
0  . 


5.1) 


The  rhombus  rules,  (R, )  and  (Rp),  which  yield  the  next  diagonal 
have  a  very  interesting  interpretation.   Recall  from  the  previous 
section  that 


o 


q 


3 


o 

In 


R 


o 


o 


From  (R-|^)  and  (R  )  and  the  rules  of  matrix  multiplication  follow  the 
result 


^k  ^  ^k-1 


o 


o 


=  ^k  +  ^k 


o   o 

q,_..e 


(R  L  ),  , 
o  o  kk 


aj^    (say) 


k+i-k    =  (^o^o^k,k+i  =  ^k  (^^y^ 


(5.2) 


It  is  easily  verified  that   RqLq   is  a  trldiagonal  matrix  of  the 


same  form  as   J   and  may  be  called  J-, 


Now  comparison  of  (5.2) 
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with  (5-1)  shows  that  the  second  diagonal  of  the  QP  scheme   (q. ,  e.) 

contains  precisely  the  elements  of  the  LR  decomposition  of  R^^L 

( =  J^ )  .   Thus 

J   Is  factorlzed  into   L  R   , 
o  o  o 

The  factors  are  multiplied  in  reverse  order  to  give   Rq^q  ^  "^1 ' 

J,   is  factorlzed  into   L.,  R^  . 
Similarly  the  third  diagonal  of  the  QD  scheme  contains  the 
elements  of  the  LR  decomposition  of   R-i  L-,  (=  Jp)  ,  and  so  on  for  all 
the  diagonals.   In  fact,  for  each  positive  integer  n  .  equations 
(5.1)  with  superscript  n  define  the  elements  a.,    p.   of  a  tri- 

diagonal  matrix   J   .   The   J   may  be  read,  as  it  were,  "between 

&  n  n 

the  (diagonal)  lines"  of  the  QD  scheme.   Clearly  J^.^.^    =  ■^n'^n^n 
and  if  the  eigenvalues   A.   of   J    (or  of  the  original  matrix  A  ) 
have  distinct  moduli  then,  as   n  — >  02 

J  — >  L    (a  triangular  matrix  with  q,  =   \^    )• 
n     00  K     K 


Thus  the  QD  scheme  implicitly  defines  a  convergent  sequence  of  tri- 
diagonal  matrices  similar  to   J 


o 


The  LR  algorithm 

This  interpretation  of  the  QD  scheme  suggests  a  similar  algo- 
rithm for  a  full  non- singular  matrix  A-^    .      By  the  LDU  theorem 
([4],  p.  20)   A-,   may  be  factorlzed  into  the  product  of  a  left 
triangular  matrix  L   and  a  right  triangular  matrix  R^  which 
has   I's   on  the  diagonal.   The  LR  transformation  is  defined  by 
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Rutishauser  ln[l4]  as  follows.   For  k  =  1,2, 
< 


Factorize   A.   Into   L,  R,  , 
k         k  k  ' 


Form  R,  L,   and  call  It  A,  , .,  . 
k  k  k-^1 

Clearly  A,  , -,  =  R,  A  R~   and  it  may  be  asked  when  this  sequence 

(A  ]   converges  to  a  triangular  matrix. 

The  answer,  see  [14],  is  contained  in  the  following 

THEOREM.   If   A^  =  U  dlag(  A^  ,  .  .  .  ,  Ajj)U"^   and 

(a)  1x^1  >  Ia^I  >  ...  >  IAj^I  , 

(b)  The  leading  principal  minors  of  U  and  U    are  non-zero, 

then  lim  A,  =  A    exists  and  is  lower  triangular  with   (A_^  ). .  =  A. 
k    00   — oo  11    1 

If  the   |a.|   are  not  distinct  then,  loosely  speaking,   A,   may 
be  said  to  converge  to  block  triangular  form.   More  precisely  suppose 
that   I  A,  I  =  |Ap|  =  ...  =  I A  I  >  U  , -,  I  •   Strictly  speaking  A^ 
may  not  exist  but,  as   k  — >  oo  ,   A,   becomes  reduced.   The  elements 
in  the  first   p  rows  may  not  converge  but  the  characteristic  poly- 
nomial of  the  leading  principal  submatrix  of  order  p   does  con- 
verge to  the  monic  polynomial  with  roots   A-,  ,  .  .  .  ,A   .   The  comple- 
mentary submatrix  has  eigenvalues  which  converge  to   A   -,,..., A^  . 
If  any  of  these  eigenvalues  have  equal  modulus  then  this  submatrix 
will  become  reduced  also.   Thus  as   k  — >  oo   the  submatrix  blocks 
which  become  isolated  along  the  diagonal  correspond  to  groups  of 
eigenvalues  of  equal  modulus.   In  the  QP  scheme,  analogously,  when 
successive   q  and  e   columns  fall  to  converge,  it  turns  out  that 
certain  polynomials,  formed  from  the   q's   and  e's   in  those 
columns  and  in  the   n  th  diagonal,  do  converge  as   n  — >  oo  . 
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In  the  important  case  of  real  matrices  with  complex  conjugate 
pairs  of  eigenvalues,   A,   may  be  expected  to  have  along  the  diagonal 
for  large   k  ,  isolated  2X  2  real  submatrices  which  yield  the 
eigenvalues  very  conveniently.   This  was  the  reason  for  stating  in 
the  previous  section  that  the  occurrence  of  eigenvalues  of  equal 
modulus  does  not  impair  the  QD  algorithm. 

Shift  of  origin  with  the  LR  transformation 

Usually  the  elements  of  A,   above  the  diagonal,  a\^.'    (i  <  j)  , 
approach  zero  like   (|a.|/|a.|)    as   k  — >  oo  .   Thus  convergence 
is  (geometric)  linear  and  depends  on  the  ratios   | A .  |  : | A  J  .   It 
happens,  quite  frequently  in  practice,  that   I^N'^''^N_iI   ^^  smaller 
than  the  other  ratios  of  consecutive  eigenvalues  and  the   (N,N) 
element  is  the  first  to  become  isolated  as   k  increases.   When 
the  remaining  elements  of  column  N  of  A,   are  zero  to  working 
accuracy  then  the   (N,N)   element  will  be  a  good  approximation  to 
A^  .   The  LR  transformation  may  then  be  applied  to  the  matrix  ob- 
tained by  omitting  row  and  column  N  .   In  this  way  the  eigenvalues 
may  be  found  in  turn,  the  order  of  the  matrix  being  reduced  at  each 
stop . 

Provided  that   s   is  considerably  closer  to  'K^     than  to  '^■^-l 

it  is  preferable  to  apply  the  transformation  to  A^-sI  rather  than 

to  A,   since  its   (N,N)   element  will  converge  to   Aj^-s   like 

(  I  A  -s  I/I  A„  -, -s  I  )    as   k  — >  oo  ;  an  improved  convergence  rate.   A 
'  N      N-1 

frequent  choice  for  s   is  the  current   (N,N)   element  but,  as 
Wilkinson  points  out  in  [22],  this  can  prevent  convergence  by  making 
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the  sequence  A,   cycle.   He  advocates  putting 

at  each  step.   This  implementation  of  the  basic  algorithm  is  very 
important,,  in  practice,  since  it  changes  the  rate  of  convergence 
from  linear  to  quadratic. 

It  might  be  objected  that  A-sI  will  be  nearly  singular  and 
hence  its  LR  factorization  will  be  liable  to  loss  of  accuracy. 
However  this  decomposition  requires  that  only  the  first  N-1   leading 
principal  minors  should  not  vanish  ([^],  p.  20).   Thus  the  trans- 
formation is  safe  numerically  even  on  A-A^I   provided  that   ^vr_-| 
is  not  too  close  to   A^  .   Yet  even  if   '^ivr_i   ai^d  A^  are  both 
small  the  algorithm  may  still  be  used.   Now  the  final  rows  and 
columns  of  the   L  and  R  factors  will  be  large  and  inaccurate 
except  for  the   (N,N)   element.   The  new  matrix  A,  .,   will  have  an 
ill-determined  last  row.   The  errors  in  this  row  will  not  impair 
the  accuracy  of  elements  in  other  rows  in  further  iterations.   When 
the   (N,N)   element  is  accepted  as   A^  -  y        s   then  this  erroneous 
last  row  is  discarded  and  the  computation  proceeds  on  the   (N-l)X(N-l) 
principal  submatrlx.   For  a  proper  discussion  of  this  same  phenomenon 
in  the  QR  algorithm  see  [6]. 

Properties  of  the  LR  transformation 

1  '5 
The  LR  decomposition  of  a  full  matrix  requires  ^  N-^  multi- 

1   3 
plications  and  the  formation  of  RL  demands  a  further  -7-  N^  .   Thus 
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the  transformation  is  not  suitable  for  full  matrices.   On  the  other 

hand  it  is  well  adapted  to  some  important  special  classes  of  matrices 

which  will  now  be  defined.   To  facilitate  the  discussion  the  writer 

proposes  the  following 

DEFINITION.   If  a.  .  =  0  for  i-j  >  m  and   j-i  >  n  then  the 

ij 

diagonal  width  of  the  matrix   [a.  .]   is   (m,n)  .   A  band  matrix  of 
width  I      is  a  matrix  whose  diagonal  width  is  {1,1)    . 

Rutishauser  [l^]  shows  that  the  width  of  band  matrices  is  pre- 
served under  the  LR  transformation  and  this  is  very  useful .   Actually 
a  more  general  result  holds. 

LEMMA .   The  diagonal  width  of  a  matrix  is  preserved  under  the 
LR  transformation. 

In  particular  the  Hessenberg  form,  diagonal  with   (N,l)   or 
(1,N)  ,  is  preserved  and  this  is  important  in  the  application  to 
general  matrices.   A  full  matrix  may  be  reduced  to  similar  Hessenberg 
form  in  a  stable  manner  with  single  precision  arithmetic  in  -^  N'^ 
multiplications,  see  Wilkinson  [l8].   The  LR  transformation  may  then 

be  applied  to  this  Hessenberg  matrix  for  which  each  iteration 

p  13 

requires   N   multiplications  as  against  -^  N"^  for  the  full  matrix. 

This  important  application  is  not  discussed  in  [l4]  but  is  exploited 

In  later  developments  of  the  LR  method. 

There  is  a  variant  of  the  basic  algorithm  given  by  Rutishauser 
for  use  with  positive  definite  symmetric  matrices  which  has  cubic 
convergence  [15].   See  also  section  8. 

In  summary  here  are  the  main  features  of  the  LR  transformation. 
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In  this  writer's  opinion,  judged  as  a  general  eigenvalue  routine. 

ADVANTAGES 

1.  The  diagonal  width  of  a  matrix  is  preserved. 

2.  The  rate  of  convergence  can  be  made  quadratic  by  judicious 
use  of  origin  shifts  (which  are  easily  coded). 

3.  Convergence  to  one  eigenvalue  (usually  the  smallest)  is 
accompanied  by  convergence  of  the  remaining  diagonal 
elements  to  the  other  eigenvalues. 

4.  The  "deflation"  technique  of  omitting  the  final  row  and 
column  when  an  eigenvalue  has  been  accepted  in  the   (N,N) 
position  is  both  safe  and  easy  to  program. 

DISADVANTAGES 

1 .  The  work  involved  in  one  iteration  on  a  full  matrix  is 
prohibitive.  However  reduction  of  the  full  matrix  to  a 
similar  Hessenberg  matrix  is  accurate  and  simple  and  then 
each  iteration  requires  less  work.  The  further  reduction 
of  a  Hessenberg  matrix  to  tridiagonal  form  is  not  without 
numerical  difficulties  (instability)  and  needs  extra  pre- 
cision in  the  computation. 

2.  The  method  is  "unstable."   More  precisely  the  LR  decompos- 
ition of  A   is  equivalent  to  performing  Gaussian  elimin- 
ation on  A  without  row  interchanges.   It  is  known  (see 
[19])  that  this  technique  can  involve  severe  loss  of 
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significant  digits  in  the  elements  of   L  and  R  .   This 
is  the  analogue  of  the  numerical  difficulties  of  the  pro- 
gressive calculation  of  the  QJ)  scheme  by  diagonals. 
One  way  of  avoiding  this  last  serious  defect  is  to  employ  a 
stable  factorization  of  a  matrix.   This  is  what  Francis  does  in  the 
QR  algorithm  [6].   Another  way  is  to  perform  the  Gaussian  elimination 
with  interchanges.   This  leads  to  a  new  algorith  of  Wilkinson  [22] 
for  real  matrices  with  real  eigenvalues.   A  variant  of  Wilkinson's 
algorithm  for  use  on  real  matrices  with  complex  eigenvalues  is  being 
studied  by  the  author.   All  three  techniques  preserve  the  Hessenberg 
form  which  is  an  important  property  for  the  practical  feasibility 
of  a  method. 
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6.   The  QR  Transformation 

For  any  matrix  A  there  exists  a  unitary  matrix  Q  and  a  right 
triangular  matrix  R  with  real,  non-negative,  diagonal  elements  such 
that 

A   =   QR  . 

Moreover  If  A   Is  non-singular  then   Q  Is  unique  and  this  Is  the 
matrix  formulation  of  the  Gram-Schmidt  orthogonallzatlon  of  the 
columns  of  A  .   Francis  In  [6]  uses  this  result  to  define  and 
develop  his  QR  transformation  as  follows:   let  A  =  A-,  ,  then  for 

k  =  1,2,... 

r 

Factorlze   A,   Into   Q,  R,   , 

^         ^  ^  (6.1) 

Form  Rif-Q;,  and  call  It  A,  ., 

*  r 

Now  A,  -  =  Q,  A  Q   and  so  the   I  A,  ]   form  a  sequence  of  unltarlly 

similar  matrices  In  contrast  to  the  LR  transformation  which  generates 

a  sequence  of  merely  similar  matrices.   It  Is  generally  believed 

that,  In  numerical  work,  unitary  transformations  are  very  stable., 

Francis  cites  Flke's  result  [5]  that  the  condition  of  an  eigenvalue 

Is  unchanged  (I.e.  not  made  worse)  when  the  matrix  undergoes  an 

orthogonal  congruence. 

The  basic  property 

Francis  proves  that  If  A   Is  non-singular  and  has  eigenvalues 
all  of  distinct  moduli  then,  as   k  — >  00  ,  the  elements  of  A, 
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below  the  principal  diagonal  tend  to  zero,  the  moduli  of  those  above 
the  diagonal  converge,  and  the  elements  on  the  principal  diagonal 
tend  to  the  eigenvalues.   If  A  has  some  eigenvalues  of  equal 
modulus  then  these  will  be  the  limits  of  the  eigenvalues  of  certain 
principal  submatrlces  of  A,   whose  elements  need  not  converge. 

In  contrast  to  LR,  the  QR  transformation  only  preserves  the 
band  width   (i,i)   of  Hermitian  matrices.   Fortunately,  however, 
the  upper  Hessenberg  form  (width   (N,l)  )  is  preserved.   Thus  a 
general  matrix  may  be  reduced  to  this  form  initially  and  the  QR 
transformation  may  then  be  applied  to  the  upper  Hessenberg  matrix 
with  great  saving  of  computation. 

In  common  with  the  LR  transformation  the  QR  transformation 
usually  yields  the  eigenvalues  ordered  in  magnitude  down  the  diagonal 
of  A    .In  practice  the  smallest  eigenvalue  usually  "appears" 
first,  located  in  the   (N,N)   element.   Thus  the  translation  tech- 
nique described  in  the  discussion  of  the  LR  method  may  be  used  here 
in  exactly  the  same  way  to  obtain  faster  convergence.   This  general- 
ized QR  transformation  is  defined  by  Francis  as  follows.   Let  A  =  A,  , 
then  for  k  =  1,2,... 

Factorize  A,  -s,  I   into   Q,,R,,   , 

k  k  k  k  (g^2) 

Form  R,  Q,  +s,  I  and  call  it  A   -  . 

The  shifts   s.   are  assumed  to  be  distinct  from,  but  approaching, 
the  smallest  eigenvalue.   There  is  some  freedom  for  skill  in  their 
choice.   Usually  one  inspects  the  roots  of  the  bottom  2X2  principal 
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submatrix.   If  they  differ  markedly  from  the  corresponding  values 
at  the  previous  Iteration  then  s,  =  0  ,   Otherwise   s.   Is  taken  as 
the  smaller  of  the  two  roots  In  magnitude. 


Note  that  A 


k+1 


C  •  •  •  ^^1^1 


I,       and  so  each  A,   Is 


similar  to  A^  .   It  would,  perhaps,  be  more  natural  to  define  ea 


ch 


A,  , ,   as   R  Q,   so  that  A,  ,,   would  be  similar  to  A.,  -  (s-, +.  .  = +s,  )I 
k+1       k  k  k+1  11      k 

The  only  objection  to  this  variant  is  that  It  would  not  allow  a 
sophisticated  trick  which  Francis  uses  on  real  matrices  and  which 
will  be  described  below.   Either  technique  may  be  used  for  complex 
matrices . 

The  formation  of  Q 

Let  A-,   be  an  upper  Hessenberg  matrix.   By  (6.2)   Q-,   Is  the 
unitary  matrix  such  that  premultlpllcatlon  of  A, -s-,1  by  Q,   elim- 
inates the  subdlagonal  elements.   There  are  only  N-1   of  these  and 

•X- 

Q-,   may  conveniently  be  taken  as  the  product  of  N-1   unitary 
matrices,  the   1  th  of  which  eliminates  element   (1+1,1)  .   Each 
may  be  taken  simply  as  the  Identity  except  for  a   2X2  principal 
submatrix  of  the  form 


6^°^  cos  e  -e^^   sin  6 


17   . 

e   sm 


15 
e   cos 


,  a+P+7+6  =  0  (mod  2tt)  ,    6   real.    (6.5) 


Such  matrices  are  the  unitary  analogues  of  Jacob! 's  two  dimensional 
rotation  matrices.   For  real  matrices  a=P=7=6=0. 

Using  these  elementary  unitary  transformations  one  step  of  the 
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QR  method  requires  approximately  4N   multiplications  as  against 

2 

N   multiplications  for  one  step  of  the  LR  method,  both  being  applied 

to  a  Hessenberg  matrix.   Stability  has  been  achieved  at  the  cost  of 
a  substantial  increase  in  the  amount  of  computation.   Of  all  the 
algorithms  considered  so  far  in  this  paper  the  QR  transformation  is 
the  first  that  can  be  used  safely  as  a  feasible  general  eigenvalue 

routine . 

The  double  QR  iteration 

Now  consider  the  important  special  case  of  real  matrices  with 
complex  eigenvalues.   In  order  to  realize  the  faster  rate  of  con- 
vergence of  the  generalized  QR  method  (using  shifts)  on  such  matrices 
complex  shifts  must  be  employed  at  some  stage.   After  that  stage  the 
matrices  will  contain  complex  elements  and  so  will  require  all  the 
extra  work  involved  in  treating  complex  matrices  (double  the  storage, 
at  least  four  times  the  work).   There  is  thus  strong  incentive  for 
seeking  an  algorithm  which  produces  a  real  matrix  from  which  the 
complex  eigenvalues  can  be  calculated  easily  as  conjugate  roots  of 
2X2   real  principal  submatrlces. 

Consider  two  steps  of  the  generalized  QR  transformation  of  a 
general  real  matrix  A^  .   From  (6.2)  it  follows  that 

A^  =  QgAgQg  =  Qg^-^Q-^Qg  (6.4) 

and  also  the  important  result  that 

Q^Q^R^Ri   =  (^-^{k^-s^DR^ 

=   (A^-SgDQ^L^l  ^^'^^ 

=      (A^-S2l)(A^-s^I)   =  G  . 
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Thus,  setting   Q  =  Q-,  Qp  and  R  =  RpR-,  ,  the  formulas  above  say 

that  A_  =  Q  A-,  Q   where,  by  uniqueness,   Q  R   Is  the   QR  factor- 
5    o  1  o       '      o  -^  o  o 

Izatlon  of  the  non-singular  matrix  G  .   It  Is  now  clear  what  should 
be  done.   Suppose  that   s-.   Is  a  complex  shift  but  that   s^   is 
forced  to  be   s-,  ,  then  G  will  be  real,   Q   will  be  its  (real) 
orthogonal  factor,   A_  will  be  real,  and  the  computation  of  complex 
A   will  be  circumvented. 

The  solution  is  therefore  to  use  the  generalized  QR  transform- 
ation with  the  extra  proviso  that  every  complex  shift  is  followed  by 
its  conjugate,  the  intermediate  complex  matrix  being  avoided  as 
indicated  above.   What  is  the  price  for  this  convenience?   Ignoring 
the  storage  of  the  auxiliary  matrix  G  it  turns  out  that  the  form- 
ation of   G  ,  its  decomposition,  and  the  formation  of  A^  {one 
double  step)  require  at  least  as  much  work  as  two  real  QR  steps. 
Apart  from  the  decided  advantage  of  keeping  the  arithmetic  real 
there  is  no  extra  gain. 

Is  there  a  better  solution?   For  a  full  real  matrix  the  answer 
is  no.   However  it  has  already  been  stated  that  it  is  a  safe  economy 
to  reduce  a  full  matrix  to  Hessenberg  form  before  using  the  QR  trans- 
formation.  Now  Francis  exploits  a  property  of  Hessenberg  matrices 

to  produce  an  algorithm  in  which  a  double  QR  step  with  complex  con- 

2  2 

jugate  shifts  requires  only   5N   multiplications  as  against   4N 

for  one  basic  real  QR  step.   This  is  most  satisfactory. 
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A  property  of  Hessenberg  matrices 

The  fact  he  uses  is  the  followingo   If  any  matrix  M  Is  unlt- 
arlly  similar  to  an  upper  Hessenberg  matrix  H  with  real,  non- 
negative,  subdlagonal  elements  then.  In  general,  the  first  column 
of  the  transformation  matrix  U  uniquely  determines  the  remaining 
columns  of  U  and  all  of  H  .   More  precisely  let 

UH  =  MU   ,  (6.6) 

then  this  equation  and  the  form  of  H  yield  a  recursive  procedure 
for  calculating  the  U.   and  the  h.  .   from  U-,   and  M  .   Subscripts 
on  matrices  denote  columns.   Let   J  =  1,2,...,N  and  at  the   j  th 
stage  take  as  known  H^,...,H._^   and  U^,...,U.  .   The  formulas 
for  H.   and  U.  ,   come  from  equating  the   J  th  column  of  each 
side  of  {6.6).      For  1  =  l,...,j   the  orthogonality  of  U^   to  the 
other  columns  of   U  gives 

h.  .   =  U?  MU.  .  (6.7) 

Using  these  values  take  all  known  quantities  to  the  right  side  of 

UH.  =  MU.   to  get 

J     J 

J 
h.  ,  .U.  T   =  MU.  -  Y~~   h,  ,U,   =  V._^T  . 

Since   U.  .   must  be  of  unit  length.  I.e.   I|U.  Jl  =  1  ,  It  follows 
J  +1  J  "^-^ 


that 


^.i,j  =  ii^-.i'i   ' 


"j.l  -    (FT^j^J.l  • 
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(6.8) 


The  procedure  comes  to  an  end  at  the  first   J  (<N)   for  which  h  .  .,  . 
vanishes.   Francis  shows  that  in  the  way  he  chooses   U-,   the  process 
will  not  end  prematurely  and  so  U-,   does  indeed  determine  the  whole 
transformation . 

Double  QR  on  a  Hessenberg  matrix 

The  application  of  this  fact  is  to  the  implication  of  (6.^)  and 

(6.5).   Take  A-,   to  be  a  real  upper  Hessenberg  matrix  (zeros  below 

the  subdiagonal)  and  recall  that  this  form  is  preserved  under  the  QR 

transformation.   So  the  orthogonal  transformation  of  A-,   into  A-, 

by  Q   is  uniquely  determined  by  the  first  column  of   Q   .   Now 

Q  R  =  G  and  R   is  right  triangular.   Therefore  column  1  of   Q 
^000^^  o 

is  just  column  1  of  G  normalized  to  unit  length.   Only  the  first 
column  of  G  is  needed  to  compute  A   and  this  column  contains  at 
most  3  non-zero  elements.   The  computation  of   G   is  avoided  thus 
saving  storage  space  and  1/12  N-^  multiplications  per  iteration. 

The  execution  of  the  method  is  the  following.   Apply  to  A-,   any 
orthogonal  congruence  the  matrix  of  which  (say   P-,  )  has  as  its  first 
column  the  normalized  first  column  of   G  .   This  will  destroy  the 
Hessenberg  form.   Now  apply  to   P-,A-,P-,   any  sequence  of  orthogonal 
congruences  (by   P_,...,P„  ,  say)  such  that 

p; ...  p;;a^p^  ...  p^  (6.9) 

is  restored  to  upper  Hessenberg  form  with  non-negative  subdiagonal 
elements.   The  one  proviso  on  the   Pp,...,P^  is  that  the  first 
column  of   P-,  .  .  .  P^  must  be  the  first  column  of   P-,  .   This  is  easily 
achieved.   Finally,  then,  if  no  subdiagonal  element  vanishes  the  new 
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•X- 

matrix  must  be  A-,  =  Q  A,  Q   since   Q   and   P.,   have  the  same  first 

3    o  1  o         o        1 

column. 


The  recursion  (6.7)   and  (6.8)  will  not  end  prematurely 

It  remains  to  show  that  no  subdiagonal  element  of  (6.9)  will 
vanish.   So  far  use  has  been  made  of  the  form  of  A^  but  not  of  the 
form  of  A   .   There  is  no  loss  of  generality  in  assuming  that  A-j^ 
has  positive  subdiagonal  elements.   (If  a  zero  occurs  then  A^   is 
reduced  and  the  two  principal  submatrlces  may  be  taken  separately.) 
Francis  proves  (Theorem  11)  that  this  assumption  ensures  that  no 
subdiagonal  element  of  A^  vanishes  provided  that   s^   is  not  an 
eigenvalue  and  that  A,   is  diagonalizable  (linear  elementary  divi- 
sors).  However  the  following  simple  argument  shows  that  this  last 
assumption  (diagonallzability )  is  not  necessary. 

Suppose  that  UH  =  A  U  where  H  is  an  upper  Hessenberg  matrix 
and  U^  =  cG-^   with   l/c  =  ||GI^||  f   0  .   Here   I^   is  the  first 
column  of   I  .   If  the  recursion  relation  (6,7)  and  (6.8)  ends 
prematurely  then  for  some   J  <  N 

J 


U.  ,   =  A,U.  -  T~  h.  .U.   =  0 
J+1      1  J   f^x   iJ  i 


Hence  there  is  a  polynomial  p(z)  ,  of  degree   j  ,  such  that 

U.^^   =  p(A^)U^   =  p(A^)GI^   =  0  . 
This  is  so  because  each  U^  ,   i  j<  j  ,  is  by  construction  the 
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vector  U,   premultlplled  by  a  polynomial  in  A   of  degree   1  . 
Now  G  Is  also  a  polynomial  In  A,   and  is  non-singular  by  hypo- 
thesis.  So   G   commutes  with   p(A,  )   and,  further,   Gp(A-,  )!-,  =  0 
implies  that   p(A,  )I   =0  .   Now  A   has  positive  subdlagonal  ele- 
ments and  so  p(A-,  )I,   has  a  non-zero   (j+1)  th  element  which  comes 
from  the  highest  degree  term  of  p  and  cannot  be  affected  by  the 
contributions  of  the  lower  degree  terms.   This  contradiction  shows 
that  this  choice  of  U   does  determine  uniquely  the  whole  trans- 
formation and  so   H  must  be  A_  . 

Francis  describes  in  detail  the  execution  of  the  algorithm  and 
the  choice  of  the   P.   for  which  he  uses  matrices  of  the  type  used 
by  Householder  to  reduce  a  symmetric  matrix  to  tridiagonal  form  [20] 
This  double  iteration  is  approximately  four  times  faster  than  a 
straightforward  computation  of  A^  via  complex  A^  . 
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7.   The  LR  Transformation  with  Interchanges 

Wilkinson  [22]  has  proposed  a  variant  of  the  LR  transformation 
which  removes  the  Instability  and  requires  no  extra  arithmetic  oper- 
ations.  This  new  technique  does  not  preserve  the  diagonal  width 
(see  section  5)  of  a  matrix  in  general  but,  fortunately,  it  does 
preserve  the  Hessenberg  form,  (diagonal  width   (N,l)   or   (1,N)  ) 
and  only  such  matrices  will  be  discussed  in  this  section.   The  method 
was  proposed  for  Hessenberg  matrices  which  are  either  real  with  real 
eigenvalues  or  complex. 

In  section  5  it  was  shown  that  one  iteration  of  the  standard  LR 
transformation  can  be  seen  as  the  result  of  premultiplylng  the 
current  matrix  A   by  a  matrix  L~    to  produce  an  upper  (or  right) 
triangular  matrix  R^  and  then  postmultlplylng  R^^  by  L^^  to 
complete  the  similarity  transformation.   L^^  Itself  is  a  product  of 
elementary  matrices  (the  identity  with  one  non-zero  subdlagonal 
element).   In  Wilkinson's  variant   L^^  is  replaced  by  a  matrix  L^^ 
which  is  not  in  general  triangular,  being  again  a  product  of  elemen- 
tary matrices  but  with  an  interchange  matrix  Inserted  between  con- 
secutive elementary  matrices.   These  Interchanges  ensure  that  no 
element  in  the  elementary  matrices  exceeds  1   in  magnitude. 

Shifts  of  origin  are  used  in  the  same  way  as  described  in  section 
5  to  increase  the  convergence  rate.   Wilkinson  reports  in  [22]  that 
this  method  is  among  the  fastest  and  most  accurate  in  his  repertoire. 
All  that  remains  is  for  a  proof  of  convergence  to  appear  in  print. 
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The  very  fact  that  seems  to  prevent  this  also  seems  to  prevent  the 
Imitation  of  the  double  Iteration  used  by  Francis  (see  section  6) 
for  real  matrices  with  complex  eigenvalues. 

How  the  Interchanges  affect  the  iteration 

Straightforward  Gaussian  elimination  used  in  solving  systems 
of  linear  equations  Is  equivalent  to  the  LR  factorization  of  the 
matrix  of  coefficients.   It  is  well  known,  see  [19]?  that  disastrous 
errors  can  be  introduced  when  the  computation  is  carried  out  with 
numbers  held  to  a  fixed  number  of  digits.   The  cure  for  this  condition 
is  to  perform  the  Gaussian  elimination  "with  Interchanges."   Before 
the  elements  below  the  main  diagonal  in  a  column  are  eliminated  two 
rows  are  interchanged  so  that  the  diagonal  element  (or  pivot)  is 
not  exceeded  in  magnitude  by  any  element  below  it.   This  method 
gives  rise  to  a  new  factorization  of  the  coefficient  matrix. 

To  see  this  in  more  detail  consider  a   5X5  upper  Hessenberg 
matrix  H,  .   Let   M.  '  denote  a   5X5  matrix  which  is  the  Identity 
except  for  a  non-zero   (1+1,1)   element.   Let   I  .,   denote  the 
identity  matrix  with  rows   a  and  b   Interchanged.   Now  the  elimin- 
ation indicated  above  may  be  written  in  matrix  notation  as 

^^4dV5c^2^2b^l^la^l   =  ^'  ' 

where  a=2orl,   b=3or2,   c=4or3,   d=5or4  accord- 
ing as  there  was  an  interchange  or  not.   When  convenient  this  equa- 
tion may  be  written  in  the  following  way 
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,-1 
where   P    =  -'-lid"'"^  "'"Ph"'"ia  '   '^®  matrices  in  parentheses,  and  hence 

,-1 
their  product   L    ,  are  left  triangular  and  so  elimination  with 

Interchanges  corresponds  to  the  factorization 

I  I  I 

H^   =   P  L  R 

in  contrast  to  the  standard  decomposition 

H-.   =   L^  R-, 
The  new  Iterate  of  H  is  given  by 

(  O)  111  I   "-'-    I   "-'-      I   I 

H^"^^   =   R  P  L    =   (L  )   (P  )   H^P  L 
in  contrast  to  the  basic  LR  iterate 

H    =   ■^1"'^1   ^   "'^1  ^1^1 

An  Important  property  of  the  basic  LR  transformation  is 
that 

H,   =  Fr-'-H^  F, 
k      k  1  k 

where   ^  S   is  the  LR  decomposition  of   (H-j^)    and  F^^  =  L^...Lj^  , 

S,  =  R  ...R-,  .   It  seems  difficult  to  find  a  usable  analogue  of  this 
k    k    1 

(2) 
property  for  the  new  algorithm.   For  example  H'     is  factored  into 

II  II  II 
P  L  R   and 
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f^)     "  "  " 

H^-'^   =   R  P  L 


L  )  ^(P  )  ^H^^^P  L 


(L  )"^(P  )"^(L  )"^(P  )  ^H-^P  L  P  L 


Also 


III'"    "I        II 


2)o'     _     /'T3't'-d'\2     _   /TT   \2 


(PLPL)(RR)   =   PLK'^'^R    =   (PLR)^   =(H^,  , 

but  the  two  matrices  on  the  left  give  neither  the  LR  nor  the  PLR 

p 
factorization  of   (H-,  )   .   This  prevents  the  straightforward  imita- 
tion of  the  double  QR  transformation,  as  is  shown  by  a  comparison 
with  the  appropriate  part  of  section  6  (formula  6.4). 

For  real  matrices  with  real  eigenvalues  It  would  suffice  to 
show  that  the  Interchanges  must  at  some  stage  become  unnecessary  for 
stability  because  after  that  stage  the  iteration  reduces  to  the 
basic  LR  transformation. 
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8.   Symmetric  Tridiagonal  Matrices 

A  popular  and  excellent  method  of  finding  the  eigenvalues  of  a 
full  real  symmetric  matrix  is  to  transform  the  matrix  into  tridiagonal 
symmetric  form  by  a  finite  sequence  of  orthogonal  transformations 
using  either  Glvens'  [7]  or  Householder's  [20]  method.   The  eigen- 
values of  the  tridiagonal  matrix  are  then  computed  successively  by 
the  Sturm  Sequence  technique  [7]  which  has  much  to  recommend  it, 
especially  if  there  are  clusters  of  close  eigenvalues.   However  there 
is  an  alternative  which  seems  to  be  from  three  to  six  times  faster 
than  the  Sturm  Sequence  algorithm  but  which  is  not  so  reliable  in 
finding  close  or  multiple  eigenvalues.   This  alternative  has  been 
noted,  among  others,  by  Wilkinson,  Ortega,  Caffrey,  and  Kaiser. 

For  a  positive  definite  real  symmetric  matrix  A-,   there  is  a 
simple  variant  of  the  LR  transformation  which  is  stable  and  preserves 

symmetry.   The  Choleskl  decomposition  of  A,   expresses  A,   as 

T 
L  L   .  where   L-,   is  a  left  triangular  matrix.   Next  Ap   is  defined 

m  T  T 

as   L,  L,   and  then  decomposed  into  L2"'^2  '  ^^  ~  "'^2''^2  ^^^   ^°  °'^* 
Rutishauser  [l4]  proved  that  under  this  transformation  A,   converges 
to  a  diagonal  matrix  as   k  — >  oo  .   This  algorithm  preserves  the 
diagonal  width  of  a  matrix  and  so  would  seem  a  sure  and  possibly 
swift  method  of  obtaining  the  eigenvalues  when  A-^^   is  tridiagonal. 
Let 
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A. 


a. 


b. 


a 


3  ^3 


3 


a 


4 


b 


^ 


a. 


From  A, 


T 


"It  -5 
it  follows  that 
,2 


L. 


"^1 

^1 

<*2 
«2 

^3 

<J4 

S. 

d. 


1 


=  ^1 


'1-1 


b./d. 


1   =  1,..  .,b,      H^ 


i,...,h 


0  , 


since  A-,   is  positive  definite  it  follows  that  the   d.  >  0  and  so 
1     ^     1 

the  square  roots   /a.-i._-,   are  certainly  real.   However  these  N 
square  roots  would  slow  down  the  numerical  execution  of  the  algorithm 
to  such  an  extent  that  it  would  not  be  competitive  with  the  Sturm 
sequence  method. 

The  modified  LR  algorithm 

The  speed  of  the  recommended  procedure  derives  from  the  fact 
that  no  square  roots  need  be  calculated.   To  see  this  consider  the 
formation  of  A    from  the  decomposition  of  A-,  .   Let 


T 


=     k. 


^1  \ 

^1   ^2   ^2 


bg  a^  5-^ 


^3  a^  b^ 
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then 

^1  =  <if  +  ^l    ,  1  =  1,...,N,   ij^  =  0 

^1   =  ^i+1-^1  '  ^   =  1,...,N-1  . 

The  Important  observation  Is  that  It  Is  not  necessary  to  know  £. 
and   d^   to  obtain  A    but  only  i.      and  dT  .   Of  course  this  means 
that   A^   Is  determined  only  up  to  the  signs  of  the   B.  .   However 
a  change  in  the  sign  of  some  (or  all)  of  the  b.   produces  a  matrix 
which  is  similar  to  Ap  via  a  diagonal  matrix  with  diagonal  elements 

±1  .   Thus  the  ambiguity  of  sign  may  be  ignored  and  A   may  be  con- 

-2 
sidered  as  determined  by  the  a.   and  b.  .   The  appropriate  modifi- 
cation of  the  straightforward  algorithm  is  Just  to  square  the  equa- 
tions giving   i.   and  b.  .   The  step  from  A.   to  A.  .,   then 
requires  only  2N  additions,   (N-1)  multiplications,  and   (N-1) 
divisions . 

Stability 

Much  has  been  made  in  this  paper  of  the  instability  of  the 
numerical  LR  transformation  for  general  matrices .   It  is  therefore 
Important  that  the  Choleskl  decomposition  of  a  positive  definite 
symmetric  matrix  is  stable.   A  thorough  discussion  of  the  matter  is 
given  by  Wilkinson  [21].   For  tridiagonal  matrices  the  equation 

0  <  d2   =  a^  -  ^,_, 

shows  that  each  \£.\      is  boiKided  by  max  a.  . 

^  i   ^ 
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It  is  also  apparent  that  the  modified  algorithm  is  applicable 
to  symmetric  matrices  which  are  not  positive  definite.   Neither  con- 
vergence to  diagonal  form  nor  stability  are  guaranteed  in  this  event 
but  for  an  Initially  positive  definite  matrix  A-,   this  modification 
does  allow  the  use  of  simple  origin  shifts  (e.g.  the   (N,N)  element) 
which  may  give  A   a  negative  eigenvalue  for  some   k  .   However  for 
positive  definite  matrices  Rutishauser  has  a  technique  [I5]  for 
choosing  the  origin  shifts  in  a  more  complicated  manner  which  yields 
cubic  convergence  to  the  diagonal  form.   It  is  not  always  advisable 
to  force  any  symmetric  matrix  to  be  positive  definite  by  use  of  a 
large  origin  shift  since  this  could  cause  the  loss  of  significant 
digits  in  the  small  elements  of  the  given  matrix. 

The  modified  QR  algorithm 

For  symmetric  trldiagonal  matrices  there  is  a  technique  whose 
stability  is  assured.   The  QR  transformation  is  greatly  simplified 
with  such  matrices  and  it  can  be  modified  so  that  no  square  roots 
need  be  taken.   Each  iteration  then  requires   5N  additions,   3N 
multiplications,  and  2N  divisions.   For  a  more  detailed  discussion 
see  Ortega  [12]. 
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